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ABSTRACT 

Curvilinear coordinate systems have been used extensively to solve 
partial differential equations on arbitrary regions. An analysis 
of truncation error in the computation of derivatives reveals why 
numerical results may be erroneous. A more accurate method of com- 
puting derivatives is presented. 


Research sponsored by NASA Langley Research Center under grant NSG 1577. 


1 . Introduction 


The use of curvilinear coordinate systems has been a major factor 
In the solution of many problems by finite difference methods. Some 
applications in the area of fluid dynamics appear in the papers by 
Godunov and Prokopov [3], Kutler [5], Steger [7], Steger and Kutler [8] 
and Thames £t. al. [9]. In this report a curvilinear coordinate system 
is a finite difference mesh with the property that each neighborhood of 
a mesh point is topologically equivalent to a rectangular mesh in the 
plane. That is, the coordinate lines can be considered as level 
curves of some one-to-one transformation. The solution of a partial 
differential equation can therefore be computed by solving the trans- 
formed equation on a rectangular region. The finite difference analog 
of the transformed equation gives a system of equations which are. to 
be solved to obtain a solution defined on the curvilinear coordinate 
system (see [9]) . 

When solving a partial differential equation on an arbitrary two- 
dimensional region, curvilinear coordinate systems can be constructed 
so that certain coordina'e lines coincide with the boundary contours 
of the region or move with the boundary in the case of a free boundary 
problem. By spacing the coordinate lines more closely in regions where 
there is a rapid change in the solution, the accuracy can often be main- 
tained with fewer mesh points than would be required with a uniform 
mesh . 

There are also difficulties which arise due to the use of a curvi- 
linear coordinate system. The transformed equation will generally be 
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more complex than the original equation and finite difference methods 
which can be used to solve the original equation may not be applicable 
to the transformed equation. Laplaces equation Is a classical example 
where efficient direct methods can be used to solve boundary value 
problems on rectangular regions. It has also been observed that the 
coordinate system can have a substantial effect on the error in the 
numerical solution of a partial differential equation. Evidence of 
this fact for one-dimensional mesh systems is amply demonstrated by 
Crowder and Dalton [2] and in a related article by Blottner and 
Roache [l]. In this report the effect of the coordinate system on 
error will be analyzed by examining the truncation error in the approx- 
imation of partial derivatives by difference expressions on the curvi- 
linear mesh. Our analysis reinforces the view that caution should be 
used in solving problems on coordinate systems with large curvature, 
rapidly changing coordinate line spacing, or with coordinate systems 
which are extremely nonorthogonal . This would be especially true 
in the case of a singular perturbation problem where the perturbation 
term might be dominated by the truncation error from other terms. 

A method of generating difference expressions is also proposed 
which incorporates some of the error terms in the traditional method 
thereby decreasing the truncation error in the derivative approxi- 
mations. Since error terms are needed, the difference expressions 
are derived from Taylor series expansions rather than by simpler method 
of differencing the transformed equations. Standard second order cen- 
tral differences have been used throughout so that an error analysis 
or an improved difference formulation could be inserted into an existing 
computer program with a minimum of modification to the code. Although 
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only quadrilateral meshes are considered In this report, a similar 
analysis could be carried out for difference formulations on triangu- 
lar meshes used by Winslow [10] and others. Related results for one- 
dimensional finite difference meshes are contained in the paper by Kelney 
de Rivas [ 4 ] and the book by Roache [6]. 

2. Difference Equations on a Curvilinear Coordinate System 

Suppose that a curvilinear coordinate system is given in a region 

3 

R of the xy-plane. Let f be a function in C (R) . Difference expres- 
sions for the first and second order partial derivatives of f can be 
obtained by transforming the curvilinear mesh to a rectangular mesh 
and applying the chain rule. The derivatives of f in terms of the 
transformed 5n-variables are related to the derivatives in the xy-plane 
by the following equations. 
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The derivatives with respect to the xy-variables can be expresced in 
terms of derivatives with respect to 5n-variables provided the Jacobian 
of the transformation does not vanish. Now consider a mesh point P with 
neighbors as indicated in Figure 1. All derivatives with respect to 
5h-variables can be approximated using classical difference operators. 
Thus we define the following expressions which replace the corresponding 
derivatives in (1) . 
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f c (P) - (f (Q) - f(R))/2 
f n (P) - (f(S) - f (T))/2 

f^(P) - f(Q) + f (R) - 2f(P) (2) 

f 5n (P) - (f(U) - f(V) + f(W) - f (X))/4 

f (P) - f(S) + f (T) - 2f(P) 

nn 

Of course these are simply the second order central differences on a 
square mesh of unit width in the £n-plane. Note that the coordinates 
x and y are also functions defined on the mesh and their derivatives 
in (1) can also be approximated by difference expressions. In this 
manner the system of equitions (1) gives rise to difference approxima- 
tions for the derivatives of f at any mesh point P. These approxima- 
tions can be used to obtain a finite difference analog of any first 
or second order partial differential equation. Unfortunately, this 
method gives no information on the accuracy of the difference approxi- 
mation. 

3. Truncation Error 

For a given curvilinear coordinate system defined in a neighbor- 
hood of P, we will associate a parameter h which will be a measure of 
the fineness of the mesh. It will be sufficient for our purposes to 
assume that the distance between any two adjacent mesh points is bounded 
above by a constant multiple of h. Thus, when f is x or y, all the 
difference expressions in (2) must be of order h, at least. 

A Taylor series expansion of the function f about the point P 


yields th^ following relations between the partial derivatives of f and 
the difference experssions in (2). 
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All remainder terms are 0(h ). In the first two equations of (3), if 
only the first order terms are retained, the remainder terms, say e^' 
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and ^ 2 * 1 wou ^ be 0(h )• In this case we would have 


3f 1/* 

3? " J <f ? y n 


W 


7 (e l' y n 


- e 


2 V 


( 4 ) 


3f 1.. 

T— = — (f X- 

3y J n 5 


£ 5 V 


J (e 2' x C 


E lV 


where 


J ■ Vn ' VC’ 


6 


The first terms on the right of (4) are the same difference approxi- 
mations for the first order partial derivatives as would be obtained 
by differencing the first two equations of (1). Consequently, the 
truncation error in the finite difference expressions for the first 
partial derivatives of f from the transformation method would be 0(h) 
provided 




- 0 ( 1 ) 


or simply 


0 ( 1 ) 


Now comparing the last three equations in (1) and (3), we see 
that the difference expressions for the second order derivatives from 
(1) would derive from the following equations which include a trunca- 
tion term. 
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Comparing (5) with (3) we note that the terms e^', and € 22 ' 

2 

would be 0(h ) since some second order terms in (3) would have to be 

combined with the remainder terms to reproduce the equations in (5) . 

3 

The determinant of the coefficient matrix is J and therefore the 
truncation error in the second derivative approximation is only 0(1) 
provided 

ItI " 0(1) • 

Jur remarks should not lead one to conclude that the equations (1) 
cannot be used to formulate accurate difference equations. The equa- 
tions obtained by this change of variable technique have been used 
extensively to solve many problems on curvilinear coordinate systems. 

It can be said that care should be used in selecting the coordinate 
system. For example, if 

“EE 1 ' |x J- ^EE 1, ly J ‘ ° <h2) ' 

then the truncation errors in the first derivative approximations 
2 

in (4) are 0(h ). No simple relation between the coefficients in the 
equation for f^ in (3) and (5) was found except for the fact that 

they would be equal if the differences were replaced by derivatives. 
Hence no general conclusions are drawn on the second derivative approx- 
imations. In the special case where the transformation is given by 
equations of the type 


x - <P(0, y ■ 'l'(n), 

then the above condition on the second order differences of the coordi- 

2 2 2 2 

nate functions implies that the truncation errors for 3 f/9x and 3 f/3y 


from (5) are 0(h). 
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Geometrically! the requirement that the second order differences 
of x and y be small means that the point P in Figure 1 and the midpoints 
of the line segments QR and ST should all be close together. This 
effectively limits the rate of change in coordinate line spacing and 
the curvature of the coordinate lines. Even for the above error esti- 
mates, it was necessary to impose limits on the quantity J. Smf.ll 

values of J also produce larger truncation errors due to the factor 1 'J 
3 

in (4) and a 1/J factor which would appear in solving the system (5) . 
Thus we see why extremely nonorthogonal coordinate systems sometimes 
give poor numerical results. 

4. Increasing Accuracy in Difference Equations 

For a general curvilinear coordinate system, more accurate dif- 
ference approximations for the partial derivatives of f can be obtained 
by using all first and second order terms of (3) . If the system (3) Is 
written in matrix form as 


A - AD + E, 

where E contains the remainder term, then the derivatives are given by 

D - A _1 A - A'^E. 


We will assume that 


det(A) 


0 ( 1 ) 


which is analogous to the previous condition on J. Since all components 

3 

of E are 0(h), the first two columns of A are 0(h) and the last three 

■> 

columns of A are 0(1 ’), it is easily shown that the first two components 
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-1 2 -1 
of A E art 0(h ) whila the last three are 0(h) . Thus A A contains 

2 

difference expressions with a truncation error of 0(h ) for the first 
order partial derivatives of f and a truncation error of 0(h) for the 
second order partial derivatives of f. 

It is difficult to determine the geometric meaning of the condition 
on det(A)» We do observe a few cases when det(A) ■ 0, however, they 
could not occur under our definition of a curvilinear coordinate system. 
?or example, if all points of Figure 1 were collnear or if P coincided 
with Q, R, S, or T, then det(A) - 0. This observation suggests that 
erroneous results could occur with coordinate systems having a value of 
det(A) which is too small. 

For many problems, the cost in computer time and storage of using 
the difference formulations of (3) rather than (A) and (5) would be 
very little. A total of 25 coefficients would be needed at each mesh 
point to compute all first and second order derivatives rather than 19 
in the latter case. The Inversion of the fifth order matrix A at each 
mesh point would not be prohibitive unless the inverse had to be recom- 
puted frequently. This would be the case if the curvilinear coordinate 
system moved during the solution of the problem as might be the case in 
solving a free boundary problem. 

5. Examples 

We now look at a couple or coordinate systems and emphasize some 

of the analysis that might be used to determine the appropriate method 

of generating difference equations. Let R be the square region give. 

by 0 <_ x,y £ 1. The coordinate lines will be parallel to the x and y 

axes and are given by the equations 

n 

x - (e N - l)/(e - 1), n ■ 0, 1, ... , N, 
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y - (e M - l)/(e - 1), m ■ 0, 1, ... , M. 


Computing the nonzero differences of x, say x^ and x^, we have the 
following approximations for large N. 


X r - 


n 

2 N 


C ' N(e - 1) 


« 


1 N 

e 


N“(e - 1) 


Similar approximations hold for the differences of y. In this example 
ve see that, for large values of M and N, the first order differences 
are much larger than the second order differen ’*3. If h is taken to 
be the minimal spacing between coordinate Hnes, then 2eh whereas 

2 2 

x^ <_ e(e - l)h . Also note that J _> 4h . A straightforward differencing 

of the transformed equations using (4) and (5) would produce difference 
approximations which are second order accurate for a first order partial 
derivative and first order accurate for a second order partial derivative. 
This is true even for the mixed derivative since (xy)^ “ x £^n an ^ 

<x %n ' ‘An ' °- 

Now a coordinate system will be defined on the same region R, but 
approximation of derivatives by (4) and (5) would not be advisable due 
to the rapid change in coordinate line spacing. Suppose the coordinate 
lines are given by 

x - (e n - 1) / (e N - 1), n - 0, 1, ... , N, 
y ■ (e ra - l)/(e^ - 1), m ■ 0, 1, ... ,M. 
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Then 


x„ ~ 


1.18 n 

N , 
e - 1 


‘55 


1.08 

N 

e - l 


n 


with similar expressions for y. In this example, the difference 
expressions for second derivatives using (5) would have a truncation 
error with terms which are nearly one forth as large as terms used in 
the difference formulation. Thus all first and second order terms 
in the series expansions (3) would be needed for accurate difference 
approximations of the derivatives of the function f. For rectangular 
coordinate systems, the determinant of the coefficient matrix A can be 
easily computed. 


. / .2, 2 1 2 W 2 12. 

det(A) - (* 5 y n ) % - ^ > <* £ - ^ > 


9 4 4 

— l^SS y nn 


> JLh 8 

-I6* 1 


The parameter r. c Id be viewed as the local minimal coordinate line 
spacing due to the large variation in coordinate line spacing through- 
out the region R. The above lower bound on det(A) is sufficient to 

guarantee that the approximation of derivatives by (3) would result 
2 

in 0(h ) and 0(h) truncation errors for derivatives of first and second 
order, respectively. 
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6. Analytic Transformations 

The analysis so far has assumed that all derivatives In the gn- 
varlables are approximated by differences. If the coordinate system 
Is generated by a simple continuous mapping function, as In the above 
example, then the derivatives of x and y can be computed and used in 
transformation equations (1) . A similar error analysis can be performed 
for the difference equations derived from this type of analytic trans- 
formation method. The previous method of generating difference equations 
will be referred to as a numerical transformation method since only the 
coordinates of the grid points are required. 

Due to the similarity in the series expansions for the analytic 
and numerical transformations, a detailed analysis of truncation error 
will not be included. It is clear from the following approximation 
formulas that the order of accuracy would depend on the orders of magni- 
tude of all derivatives of x and y with respect to g and n- 

For the first order central differences in the computational region, 
we will assume that f, as a function of g and n, is sufficiently smooth 
so that 


f 


£ 


3f 1 jiff 

95 6 S£ 3 . 


( 6 ) 


The right hand sides will now be replaced by derivatives with respect to 
th 2 physical variables. In order to have a comparison with numerical 
transformations, only the first and second order derivatives are retained. 
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The Inclusion of additional terms in (6) would introduce additional 
higher order terms in the coefficients of the derivatives of f. 
Approximations for the second order differences are derived in a 
similar fashion. If the assumption is made that 

, . ifj 4.1. ifl f _ 3 2 f 1 3-f 

K H 2 12 35* '*> >«n - 6^ + 

then the following estimates hold after dropping all but the first and 
second order terms. 
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Again we observe that the traditional difference formulations for the 
derivatives of f are accurate only if the higher order derivatives of 
the coordinate functions x and y become progressively smaller. 

In the above error analysis for both numerical and analytic trans- 
formations, it was assumed that certain remainder terms in the Taylor 
series expansions were negligible. This implies the third order deriva- 
tives of f with respect to the physical variables must be bounded and 
these bounds depend on the coordinate line spacing. A comparison of 
the above approximations also suggest that the numerical computation of 
the derivative coefficients would be preferred due to the appearance of 
error in the first derivative terms when the coefficients are computed 
analytically. On the other hand, for the numerical transformation, the 

g2f 

coefficients of the second derivative terms in the expression for - ■ ■ ■ 
appear to have little relation to the coefficients encountered in the 
usual change of variables formulation. The differential analogs of the 
corresponding difference expressions for the coefficients would, however, 
be equal. 

Equivalent estimates would hold for central differences for functions 
of three variables. These expressions are omitted since they are even 
more length than those given above. 
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